Chapter 5 Community composition

load("resources/data.Rdata")

5.1 Taxonomy overview

5.1.1 Stacked barplot

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  mutate(individual=factor(individual,levels=c("Sg1","Sg2","Sg3","Sg4","Sg5","Sg6","Sg7","Sg8","Sg9","Sg10"))) %>%
  ggplot(., aes(x=sample,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(. ~ individual,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

5.1.2 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun),sd=sd(relabun)) %>%
    arrange(-mean) %>%
    tt()
tinytable_hchg3d8rycd06ridhf4s
phylum mean sd
p__Pseudomonadota 0.4043587138 0.464510230
p__Bacillota_A 0.2355141069 0.244605492
p__Bacteroidota 0.2175456020 0.236021607
p__Campylobacterota 0.0994533579 0.304943152
p__Bacillota 0.0188739477 0.033201006
p__Desulfobacterota 0.0116207168 0.017287755
p__Verrucomicrobiota 0.0059970322 0.008688895
p__Bacillota_B 0.0035458974 0.008344688
p__Bacillota_C 0.0025972812 0.004149428
p__Cyanobacteriota 0.0004933441 0.002206302 
phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    select(phylum) %>%
    pull()

phylum_summary %>%
    filter(phylum %in% phylum_arrange) %>%
    mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
    ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
        scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        geom_jitter(alpha=0.5) + 
        theme_minimal()

## Taxonomy boxplot

5.1.3 Family

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_summary %>%
    group_by(family) %>%
    summarise(mean=mean(relabun),sd=sd(relabun)) %>%
    arrange(-mean) %>%
    tt()
tinytable_hkucymqic0jy7jgmdegz
family mean sd
f__Enterobacteriaceae 0.4037464761 0.465021606
f__Lachnospiraceae 0.2133481440 0.225431409
f__Bacteroidaceae 0.1153186642 0.123382341
f__Helicobacteraceae 0.0994533579 0.304943152
f__Rikenellaceae 0.0486474961 0.064078311
f__Tannerellaceae 0.0378215379 0.044036446
f__Marinifilaceae 0.0157579038 0.021580609
f__Desulfovibrionaceae 0.0116207168 0.017287755
f__Coprobacillaceae 0.0105137444 0.028390737
f__Oscillospiraceae 0.0085479249 0.014371104
f__Erysipelotrichaceae 0.0072817637 0.012422508
f__Akkermansiaceae 0.0059970322 0.008688895
f__Ruminococcaceae 0.0039836005 0.013192296
f__Acutalibacteraceae 0.0038692915 0.011317322
f__ 0.0036487884 0.005954936
f__Peptococcaceae 0.0035458974 0.008344688
f__UBA3700 0.0020129408 0.004923646
f__CAG-508 0.0015982064 0.002962023
f__Anaerovoracaceae 0.0011024918 0.002537585
f__UBA660 0.0010784396 0.002763385
f__CAG-239 0.0006122376 0.001950787
f__Gastranaerophilaceae 0.0004933441 0.002206302 
family_arrange <- family_summary %>%
    group_by(family) %>%
    summarise(mean=sum(relabun)) %>%
    arrange(-mean) %>%
    select(family) %>%
    pull()

family_summary %>%
    left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==sample)) %>%
    filter(family %in% family_arrange[1:20]) %>%
    mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        #geom_boxplot() +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal()

5.1.4 Genus

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_summary %>%
    group_by(genus) %>%
    summarise(mean=mean(relabun),sd=sd(relabun)) %>%
    arrange(-mean) %>%
    tt()
tinytable_a8199pbhgbfovywmbrtf
genus mean sd
g__Hafnia 0.2932325053 0.436327554
g__Salmonella 0.1042947262 0.306889524
g__Bacteroides 0.0679517227 0.073152380
g__Alistipes 0.0470015324 0.060867750
g__Phocaeicola 0.0401804345 0.047455237
g__Parabacteroides 0.0350520755 0.041334502
g__Hungatella 0.0270817008 0.043700922
g__Roseburia 0.0254184507 0.065658830
g__Clostridium_Q 0.0241488625 0.034458097
g__Ventrimonas 0.0214487733 0.028443966
g__Enterocloster 0.0200934974 0.030265361
g__Odoribacter 0.0157579038 0.021580609
g__Acetatifactor 0.0148914484 0.031889132
g__Eisenbergiella 0.0147363777 0.021995069
g__Bilophila 0.0098378050 0.015392491
g__Thomasclavelia 0.0076567822 0.026812811
g__Hungatella_A 0.0062116286 0.014527410
g__OM05-12 0.0060236930 0.008143623
g__Akkermansia 0.0059970322 0.008688895
g__Ventrisoma 0.0057439870 0.018833487
g__Breznakia 0.0047686166 0.010282743
g__Kineothrix 0.0047578797 0.014794492
g__Escherichia 0.0039316851 0.014505764
g__JAAYNV01 0.0033683384 0.010551754
g__Anaerotruncus 0.0029580653 0.013228870
g__Parabacteroides_B 0.0027694623 0.005044226
g__JAAWBF01 0.0026779996 0.005456383
g__UBA7185 0.0025995472 0.007619166
g__Dielma 0.0025131471 0.003740659
g__14-2 0.0024519108 0.005213359
g__Klebsiella 0.0022875595 0.006505467
g__Velocimicrobium 0.0022666655 0.004653853
g__Pseudoflavonifractor 0.0021662766 0.006472208
g__Lawsonibacter 0.0019945045 0.002717223
g__Desulfovibrio 0.0017829118 0.003516870
g__NSJ-51 0.0017091442 0.003617896
g__Tidjanibacter 0.0016459637 0.005192887
g__Merdicola 0.0015982064 0.002962023
g__Lachnoclostridium_A 0.0014855570 0.003070051
g__Lacrimispora 0.0014491692 0.003818055
g__Murimonas 0.0012522702 0.003169218
g__Onthovicinus 0.0011673635 0.005220608
g__Bacteroides_G 0.0011628140 0.003033977
g__Emergencia 0.0011024918 0.002537585
g__Faecimonas 0.0010784396 0.002763385
g__MGBC101980 0.0010641277 0.003493127
g__Novisyntrophococcus 0.0010323253 0.002311560
g__CAZU01 0.0006122376 0.001950787
g__Limenecus 0.0004933441 0.002206302 
genus_arrange <- genus_summary %>%
    group_by(genus) %>%
    summarise(mean=sum(relabun)) %>%
    filter(genus != "g__")%>%
    arrange(-mean) %>%
    select(genus) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    pull()

genus_summary %>%
    left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
    left_join(sample_metadata,by=join_by(sample==sample)) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    filter(genus %in% genus_arrange[1:20]) %>%
    mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
        scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
        #geom_boxplot() +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal()

5.2 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Sample type",
          breaks=c("cloaca","feces"),
          labels=c("Cloaca","Faeces"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="Sample type",
          breaks=c("cloaca","feces"),
          labels=c("Cloaca","Faeces"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric, scales = "free", ncol=4) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

alpha_div %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  lmerTest::lmer(richness ~ type + (1 | individual), data = ., REML = FALSE) %>%
  broom.mixed::tidy() %>%
  tt()
boundary (singular) fit: see help('isSingular')
tinytable_g28n7vz8kq04tom1n6sg
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 1.600000 2.987725 0.5355245 20 5.981925e-01
fixed NA typefeces 61.300000 4.225281 14.5079106 20 4.444567e-12
ran_pars individual sd__(Intercept) 0.000000 NA NA NA NA
ran_pars Residual sd__Observation 9.448016 NA NA NA NA 
alpha_div %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  lmerTest::lmer(neutral ~ type + (1 | individual), data = ., REML = FALSE) %>%
  broom.mixed::tidy() %>%
  tt()
boundary (singular) fit: see help('isSingular')
tinytable_zhay84pyo40vnt5jxcy3
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 1.300518e+00 1.359874 0.9563518 20 3.503127e-01
fixed NA typefeces 3.294989e+01 1.923152 17.1332730 20 2.024047e-13
ran_pars individual sd__(Intercept) 4.971903e-09 NA NA NA NA
ran_pars Residual sd__Observation 4.300298e+00 NA NA NA NA 
alpha_div %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  lmerTest::lmer(phylogenetic ~ type + (1 | individual), data = ., REML = FALSE) %>%
  broom.mixed::tidy() %>%
  tt()
boundary (singular) fit: see help('isSingular')
tinytable_hf86gfsu98xygfwhz7zx
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 1.2384408 0.2335543 5.302583 20 3.441968e-05
fixed NA typefeces 3.1026737 0.3302956 9.393626 20 8.961863e-09
ran_pars individual sd__(Intercept) 0.0000000 NA NA NA NA
ran_pars Residual sd__Observation 0.7385635 NA NA NA NA 
alpha_div %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  lmerTest::lmer(functional ~ type + (1 | individual), data = ., REML = FALSE) %>%
  broom.mixed::tidy() %>%
  tt()
boundary (singular) fit: see help('isSingular')
tinytable_ltktuosf3gevro2iqi0d
effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 1.052834e+00 0.03546038 29.690432 20 5.139978e-18
fixed NA typefeces 4.136317e-01 0.05014855 8.248128 20 7.248330e-08
ran_pars individual sd__(Intercept) 3.057788e-11 NA NA NA NA
ran_pars Residual sd__Observation 1.121356e-01 NA NA NA NA

5.3 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)
#Richness
betadisper(beta_q0n$C, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.07339 0.073391 0.6442    999  0.428
Residuals 18 2.05080 0.113933                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        cloaca feces
cloaca         0.423
feces  0.43268      
adonis2(beta_q0n$C ~ type, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
tinytable_xrswy83to1xy3wl6fbgi
term df SumOfSqs R2 statistic p.value
type 1 3.005780 0.4664694 15.73752 0.002
Residual 18 3.437902 0.5335306 NA NA
Total 19 6.443682 1.0000000 NA NA 
#Neutral diversity
betadisper(beta_q1n$C, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01774 0.01774 0.1446    999  0.694
Residuals 18 2.20802 0.12267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        cloaca feces
cloaca         0.703
feces  0.70818      
adonis2(beta_q1n$C ~ type, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
tinytable_82hj1v7r9oh04w7mfsjk
term df SumOfSqs R2 statistic p.value
type 1 2.322537 0.362287 10.22586 0.002
Residual 18 4.088228 0.637713 NA NA
Total 19 6.410765 1.000000 NA NA 
#Phylogenetic diversity
betadisper(beta_q1p$C, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.07195 0.071945 1.0792    999  0.358
Residuals 18 1.20001 0.066667                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        cloaca feces
cloaca         0.449
feces  0.31264      
adonis2(beta_q1p$C ~ type, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
tinytable_ai77kq006fwiybxz0e3i
term df SumOfSqs R2 statistic p.value
type 1 2.893508 0.6789794 38.07116 0.002
Residual 18 1.368047 0.3210206 NA NA
Total 19 4.261555 1.0000000 NA NA 
#Functional diversity
betadisper(beta_q1f$C, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.15953 0.159526 7.2962    999  0.015 *
Residuals 18 0.39355 0.021864                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         cloaca feces
cloaca          0.013
feces  0.014616      
adonis2(beta_q1f$C ~ type, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
tinytable_3srr1rsuusadixf79dza
term df SumOfSqs R2 statistic p.value
type 1 3.1347963 0.8009941 72.44958 0.002
Residual 18 0.7788359 0.1990059 NA NA
Total 19 3.9136321 1.0000000 NA NA

5.3.1 Neutral diversity plot

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, verbosity = FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  mutate(individual=factor(individual, levels=c("Sg1","Sg2","Sg3","Sg4","Sg5","Sg6","Sg7","Sg8","Sg9","Sg10"))) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = type, shape = individual)) +
    scale_color_manual(name="Sample type",
          breaks=c("cloaca","feces"),
          labels=c("Cloaca","Faeces"),
          values=c("#e5bd5b", "#6b7398")) +
    scale_shape_manual(values = 1:10) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(shape="Individual")

5.3.2 Functional diversity plot

beta_q1f$C %>%
  vegan::metaMDS(., trymax = 500, k = 2, verbosity = FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(individual=factor(individual, levels=c("Sg1","Sg2","Sg3","Sg4","Sg5","Sg6","Sg7","Sg8","Sg9","Sg10"))) %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = type)) +
    scale_color_manual(name="Sample type",
          breaks=c("cloaca","feces"),
          labels=c("Cloaca","Faeces"),
          values=c("#e5bd5b", "#6b7398")) +
    scale_shape_manual(values = 1:10) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )